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STABILIZATION IN RELATION TO WAVENUMBER IN HDG 

METHODS 


J. GOPALAKRISHNAN, S. LANTERI, N. OLIVARES, AND R. PERRUSSEL 

Abstract. Simulation of wave propagation through complex media relies on proper 
understanding of the properties of numerical methods when the wavenumber is real 
and complex. Numerical methods of the Hybrid Discontinuous Galerkin (HDG) type 
are considered for simulating waves that satisfy the Helmholtz and Maxwell equations. 
It is shown that these methods, when wrongly used, give rise to singular systems for 
complex wavenumbers. A sufficient condition on the HDG stabilization parameter for 
guaranteeing unique solvability of the numerical HDG system, both for Helmholtz and 
Maxwell systems, is obtained for complex wavenumbers. For real wavenumbers, results 
from a dispersion analysis are presented. An asymptotic expansion of the dispersion 
relation, as the number of mesh elements per wave increase, reveal that some choices of 
the stabilization parameter are better than others. To summarize the findings, there 
are values of the HDG stabilization parameter that will cause the HDG method to 
fail for complex wavenumbers. However, this failure is remedied if the real part of the 
stabilization parameter has the opposite sign of the imaginary part of the wavenumber. 
When the wavenumber is real, values of the stabilization parameter that asymptoti¬ 
cally minimize the HDG wavenumber errors are found on the imaginary axis. Finally, 
a dispersion analysis of the mixed hybrid Raviart-Thomas method showed that its 
wavenumber errors are an order smaller than those of the HDG method. 


Background 

Wave propagation through complex structures, composed of both propagating and 
absorbing media, are routinely simulated using numerical methods. Among the vari¬ 
ous numerical methods used, the Hybrid Discontinuous Galerkin (HDG) method has 
emerged as an attractive choice for such simulations. The easy passage to high order 
using interface unknowns, condensation of all interior variables, availability of error es¬ 
timators and adaptive algorithms, are some of the reasons for the adoption of HDG 
methods. 

It is important to design numerical methods that remain stable as the wavenumber 
varies in the complex plane. For example, in applications like computational lithography, 
one finds absorbing materials with complex refractive index in parts of the domain of 

Key words and phrases. HDG, Raviart-Thomas, dispersion, dissipation, absorbing material, com¬ 
plex, wave speed, optimal, stabilization, Helmholtz, Maxwell, unisolvency. 

JG and NO were supported in part by the NSF grant DMS-1318916 and the AFOSR grant FA9550- 
12-1-0484. NO gratefully acknowledges support in the form of an INRIA internship where discussions 
leading to this work originated. All authors wish to thank INRIA Sophia Antipolis Mediterranee for 
hosting the authors there and facilitating this research. 

1 



2 


J. GOPALAKRISHNAN, S. LANTERI, N. OLIVARES, AND R. PERRUSSEL 


simulation. Other examples are furnished by meta-materials. A separate and important 
reason for requiring such stability emerges in the computation of resonances by iterative 
searches in the complex plane. It is common for such iterative algorithms to solve a 
source problem with a complex wavenumber as its current iterate. Within such algo¬ 
rithms, if the HDG method is used for discretizing the source problem, it is imperative 
that the method remains stable for all complex wavenumbers. 

One focus of this study is on complex wavenumber cases in acoustics and electromag¬ 
netics, motivated by the above-mentioned examples. Ever since the invention of the 
HDG method in [4], it has been further developed and extended to other problems in 
many works (so many so that it is now impractical to list all references on the subject 
here). Of particular interest to us are works that applied HDG ideas to wave propaga¬ 
tion problems such as [6, 8, 10, 11, 12, 13, 14], We will make detailed comparisons with 
some of these works in a later section. However, none of these references address the 
stability issues for complex wavenumber cases. While the choice of the HDG stabiliza¬ 
tion parameter in the real wave number case can be safely modeled after the well-known 
choices for elliptic problems [5], the complex wave number case is essentially different. 
This will be clear right away from a few elementary calculations in the next section, 
which show that the standard prescriptions of stabilization parameters are not always 
appropriate for the complex wave number case. This then raises further questions on 
how the HDG stabilization parameter should be chosen in relation to the wavenumber, 
which are addressed in later sections. 

Another focus of this study is on the difference in speeds of the computed and the 
exact wave, in the case of real wavenumbers. By means of a dispersion analysis, one can 
compute the discrete wavenumber of a wave-like solution computed by the HDG method, 
for any given exact wavenumber. An extensive bibliography on dispersion analyses 
for the standard finite element method can be obtained from [1, 7]. For nonstandard 
finite element methods however, dispersion analysis is not so common [9], and for the 
HDG method, it does not yet exist. We will show that useful insights into the HDG 
method can be obtained by a dispersion analysis. In multiple dimensions, the discrete 
wavenumber depends on the propagation angle. Analytic computation of the dispersion 
relation is feasible in the lowest order case. We are thus able to study the influence 
of the stabilization parameter on the discrete wavenumber and offer recommendations 
on choosing good stabilization parameters. The optimal stabilization parameter values 
are found not to depend on the wavenumber. In the higher order case, since analytic 
calculations pose difficulties, we conduct a dispersion analysis numerically. 

We begin, in the next section, by describing the HDG methods. We set the stage for 
this study by showing that the commonly chosen HDG stabilization parameter values 
for elliptic problems are not appropriate for all complex wavenumbers. In the subse¬ 
quent section, we discover a constraint on the stabilization parameter, dependent on the 
wavenumber, that guarantees unique solvability of both the global and the local HDG 
problems. Afterward, we perform a dispersion analysis for both the HDG method and 
a mixed method and discuss the results. 
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Methods of the HDG type 


We borrow the basic methodology for constructing HDG methods from [4] and apply it 
to the time-harmonic Helmholtz and Maxwell equations (written as first order systems). 
While doing so, we set up the notations used throughout, compare the formulation 
we use with other existing works, and show that for complex wavenumbers there are 
stabilization parameters that will cause the HDG method to fail. 


Undesirable HDG stabilization parameters for the Helmholtz system. We 

begin by considering the lowest order HDG system for Helmholtz equation. Let A; be a 
complex number. Consider the Helmholtz system on fl c M 2 with homogeneous Dirichlet 
boundary conditions, 


ikU + V<P = 0, 

in H, 

(la) 

ik<P + SJ-U = /, 

in H, 

(lb) 

<k> = 0, 

on dkl, 

(lc) 


where / e L 2 {Vt). Let Th denote a square or triangular mesh of disjoint elements K , so 
fl = u KeTh K, and let Th denote the collection of edges. The HDG method produces an 
approximation (-5,4,4) to the exact solution where & denotes the trace of <P 

on the collection of element boundaries <974 The HDG solution (5,4, 4) is in the finite 
dimensional space 14 x H4 x M h defined by 

14 = {h e (L 2 (H)) 2 : v\ K €V(K), \/K € %,} 

W h = {^eL 2 (n) : 4| K eW(K), VW € 77,} 

M h = {4 e L 2 ( U F) : 41 f e M(F), VF e T h and 4|an = 0}, 

F*T h 

with polynomial spaces V(K), W(K), and M(F ) specihed differently depending on 
element type: 


Triangles 
V(K) = (V P (K)) 
W(K)=V P (K) 
M(F)=V P (F) 


Squares 

V(K) = (Q P (K)) 
W(K) = Q P (K) 
M(F)=V P (F). 


Here, for a given domain D, V P (D) denotes polynomials of degree at most p, and Q P (D ) 
denotes polynomials of degree at most p in each variable. 

The HDG solution satisfies 


E 

KeT h 

ik(u, v) K - (4, V- v ) K + {4>,v- n)di< = 0, 

(2a) 

E 

KeT h 

-(v-u, 4) A ' + (r4,4W - (t 4,4W -?fc(4,4)tf = -(/,4)n> 

(2b) 

E 

KeT h 

(u-n + r{4>-4>)A)dK = 0, 

(2c) 
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Reference 

Their notations and equations 

Connection to our formulation 

[6] 

Helmholtz case 

%i + % = 8 

'W,r ,A '.;r (l 

■'ter 

T i°i = kT 
iku [6] = 0 

%] = 

[10] 

Helmholtz case 

ikq m + V« [10] - 0 

*WU[io] + V-g[io] = ° 

^[101 ’ ™ = Vi ’ ™ + WVi _ “[ 101 ^ 

r [io] = r 

M [10] = ^ 

V] = 

[12] 

2D Maxwell case 

^[i 2 ]^[i 2] - Vx ^ r [i 2] = ° 

% 12 frH [12] + VxE [l2] = 0 
^ [12] ~ ^[12] + r [12]^[12] ~ ^[12]^ 

j~Sr 

r i 12 l = V ^r T 

w [12] = UJ \/ £ 0f i 0 

~ ^ r E ’ E W ■ ^ TT r E 

[14] 

Maxwell case 

^[14]-^ X V] = 5 

Vxu5 [14] - £ u; 2 u [1 4] = 0 

iS [ 14] = >4] + T [1 4 ] ( V] _ V] )X ” 

„ jeoj 2 

Vr*V h r 

fiw = -ikH , with k = (jy/JH, 

Vi = E 


Table 1. Comparison with some HDG formulations in other papers. No¬ 
tations in the indicated external references are used after subscripting 
them by the reference number. Notations without subscripts are those 
defined in this paper. 


for all v e 14, 0 e R0, and 0 € M h . The last equation enforces the conservativity of the 
numerical flux 


u ■ n - u ■ n + r(0 - 0). (3) 

The stabilization parameter r is assumed to be constant on each dK. We are interested 
in how the choice of r in relation to k affects the method, especially when k is complex 
valued. Comparisons of this formulation with other HDG formulations for Helmholtz 
equations in the literature are summarized in Table 1. 

One of the main reasons to use an HDG method is that all interior unknowns (u, 0) 
can be eliminated to get a global system for solely the interface unknowns (0). This is 
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possible whenever the local system 

ik{u,v) K - {(j), V-v) K = 
-(W-u,^) K - {r(j), ip) eK - ip) K = 


-{4>,v-n) dK , 
-{r&'ipjdK 
~ (/Y)x, 


Vu € V(K), 

€ W(K ), 


(4a) 

(4b) 


is uniquely solvable. (For details on this elimination and other perspectives on HDG 
methods, see [4].) In the lowest order (p = 0) case, on a square element K of side length 
h, if we use a basis in the following order 



Y 


'o' 

u x = 

0 

, U 2 = 

1 


then the element matrix for the system (4) is 


1, on K, 


M 


'ik h 2 0 0 

ik h? 0 

0 0-4 hr - ik h 2 


This shows that if 

4t = - ikh , (5) 

then M is singular, and so the HDG method will fail. The usual recipe of choosing t = 1 
is therefore inappropriate when k is complex valued. 


Intermediate case of the 2D Maxwell system. It is an interesting exercise to 
consider the 2D Maxwell system before going to the full 3D case. In fact, the HDG 
method for the 2D Maxwell system can be determined from the HDG method for the 
2D Helmholtz system. The 2D Maxwell system is 

ik£ - Vx % = - J, 
ikH + Vx £ - 0, 

where J e L 2 (£l), and the scalar curl Vx • and the vector curl Vx • are defined by 

Vx n = ft U 2 - d 2 U x - V- R(H), V*£= (d 2 £, -di£) = R(y£). 

Here R(vi,v 2 ) = (v 2 ,-vi) is the operator that rotates vectors counterclockwise by +n/2 
in the plane. Clearly, if we set f = -Rfi, then (6) becomes 

ik£ + V- r = - J, 

-ikr + RR\j£ = 0, 

which, since RRv = -v (rotation by 7r), coincides with (1) with <P - £, U - f, and / = - J. 
This also shows that the HDG method for Helmholtz equation should yield an HDG 
method for the 2D Maxwell system. We thus conclude that there exist stabilization 
parameters that will cause the HDG system for 2D Maxwell system to fail. 


(6a) 

(6b) 
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To examine this 2D HDG method, if we let H and E denote the HDG approximations 
for Rf and £, respectively, then the HDG system (2) with u and (j) replaced by -RH 
and E, respectively, gives 

Y ~( E , Vx w) K + (E, n x w) dK - ik(H, w) K = 0, 

KeT h 

Y ik(E, ip) K - (Vx H 7 ij)) K + (t(E-E), 'll ))ok = 

KzT h 

Y ( R H ■ n, )dK = o, 

KsT h 

for all w e R{V h )^ e Wh and ^ € M h . We have used the fact that -(RH) ■ n = H -t, 
where t - Rn the tangent vector, and we have used the 2D cross product defined by 
v x fi = v ■ t. In particular, the numerical flux prescription (3) implies 

-RH ■ n = -RH ■ n + r(E - E), 

where RH denotes the numerical trace of RH. We rewrite this in terms of H and E, to 
obtain 

H - t = H -t + r(E - E). 

One may rewrite this again, as 

Hxn = Hxn + r(E-E). (7) 

This expression is notable because it will help us consistently transition the numerical 
flux prescription from the Helmholtz to the full 3D Maxwell case discussed next. A 
comparison of this formula with those in the existing literature is included in Table 1. 

The 3D Maxwell system. Consider the 3D Maxwell system on ficR 3 with a perfect 
electrically conducting boundary condition, 


lk£ - Vx R - - J, 

in 0, 

(8a) 

ikH + Vx S - 0, 

in D, 

(8b) 

fix E - 0, 

on d Q, 

(8c) 


where J e (L 2 (Q)) 3 . For this problem, Th denotes a cubic or tetrahedral mesh, and Eh 
denotes the collection of mesh faces. The HDG method approximates the exact solution 
(£. 77, £) by the discrete solution (E, H, E) € W x W x W• The discrete spaces are defined 
by 

Y h = {ve(L 2 (n )) 3 : v\ K € Y(K), VA e T h ) 

N h = {fj e (L 2 (E h )) 3 : fj\ F e N(F), VF e T h and f)\dn = 6}, 
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with polynomial spaces Y(K ) and N(F ) specified by: 

Tetrahedra Cubes 

Y(K) = (V P (K)) 3 Y(K) = (Q P (K)) 3 

N(F) = {fj € (V P (F)) 3 : fj -n = 0} N(F) = {fj e (Q P (F)) 3 : ? 7 -fi = 0} 

Our HDG method for (8) is 


E 

KeT h 

ik(E, v) K ~ ( Vx H,v) k + ((H-H) x n, v) dK = -(J,fi)o, 

Vfi € Y h , 

E 

KeT h 

-(A, Vx fi))x + (A 1 , fi x ~ ik(H,w)K = 0, 

Vfi) € Y hi 

E 

KeT h 

(ifi x n, w) dK = 0 , 

VD e A^ h , 

where, in 

analogy with (7), we now set numerical flux by 



H x n - H x n + t(E - E) t , 

(9) 

where (E 

- E) t denotes the tangential component, or equivalently 



Hxn = Hxn + r{n x (E - E)) x n. 


Note that the 2D system ( 6 ) is obtained from the 3D Maxwell system ( 8 ) by assuming 


symmetry in 23 -direction. Hence, for consistency between 2D and 3D formulations, we 
should have the same form for the numerical flux prescriptions in 2D and 3D. 

The HDG method is then equivalently written as 

Y ik(E,v) K - (Vx H,v) k + (r(E-E) x n, fi x n) dK = (10a) 

I<cT h 

^ ~(E, Vx w)k + (E, n x w)g K - ik(H, w)k - 0, (10b) 

KeTh 

^ (H + rn x (E - E), w x n)dK = 0, (10c) 

KeT h 

for all v,w £ Y h , and w e Nh- For comparison with other existing formulations, see 
Table 1. 

Again, let us look at the solvability of the local element problem 

ik(E,v)K ~ (Vx H,v)k + (tE xn, fix u)qk = ( tE x fit, fi x n)g K 

~(J,v) K , (Ha) 

-(F,Vxi u) K -ik{H,w) K = -( E,n x fi5)ax, (lib) 

for all fi,rfi € Y(K). In the lowest order (p = 0) case, on a cube element K of side length 
h, if we use a basis in the following order 
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then the 6 x 6 element matrix for the system ( 11 ) is 

(4 h 2 r + ikh 3 )I 3 0 

0 ~{ikh 3 )h j * 

where denotes the 3x3 identity matrix. Again, exactly as in the Helmholtz case - 
cf. (5) - we find that if 

4r = -ikh, (13) 

then the local static condensation required in the HDG method will fail in the Maxwell 
case also. 

Behavior on tetrahedral meshes. For the lowest order (p = 0) case on a tetrahe¬ 
dral element, just as for the cube element described above, there are bad stabilization 
parameter values. Consider, for example, the tetrahedral element of size h defined by 

K - {x e M 3 : Xj ^ 0 Vj, xi + x 2 + £3 5 $ h}, (14) 

with a basis ordered as in (12). The element matrix for the system (11) is then 

(2^3 + 6 )h 2 r + ikh 3 h 2 r -\f?>h 2 T 0 0 0 

-\/3 h 2 r (2\/3 + 6 )h 2 r + ikh 3 -\/Zh 2 T 0 0 0 

-y/3 h 2 r -\/3h 2 T Ah 2 r + ikh 3 0 0 0 

0 0 0 -ikh 3 0 0 

0 0 0 0 -ikh 3 0 

0 0 0 0 0 -ikh 3 _ 

We immediately see that the rows become linearly dependent if 

(3v^3 + 6 )r = -ikh. 

Hence, for r = -ikh/( 3\/3 + 6 ), the HDG method will fail on tetrahedral meshes. 

For orders p ^ 1, the element matrices are too complex to find bad parameter values 
so simply. Instead, we experiment numerically. Setting r = -i, which is equivalent to 
the choice made in [14] (see Table 1), we compute the smallest singular value of the 
element matrix M (the matrix of the left hand side of (11) with K set by (14)) for a 
range of normalized wavenumbers kh. Figures la and lb show that, for orders p = 1 and 
p = 2 , there are values of kh for which r = -i results in a singular value very close to 
zero. Taking a closer look at the first nonzero local minimum in Figure la, we find that 
the local matrix corresponding to normalized wavenumber kh » 7.49 has an estimated 
condition number exceeding 3.9 x 10 15 , i.e., for all practical purposes, the element matrix 
is singular. To illustrate how a different choice of stabilization parameter r can affect 
the conditioning of the element matrix, Figures lc and Id show the smallest singular 
values for the same range of kh, but with r = 1. Clearly the latter choice of r is better 
than the former. 

From another perspective, Figure le shows the smallest singular value of the element 
matrix as r is varied in the complex plane, while fixing kh to 1. Figure If is similar 
except that we fixed kh to the value discussed above, approximately 7.49. In both cases, 
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(A) p= 1, r = -i 



(c) p= 1, r= 1 


(b) p = 2, t - -i 



(d) p = 2, t = 1 



(e) kh = 1, p = 1 


(f) fc/i ^ 7.49, p = 1 



(g) kh = 1 + z, p = 1 


(h) kh pa 7.49(1 + 1), p = 1 


Figure 1. The smallest singular values of a tetrahedral HDG element matrix 




















10 


J. GOPALAKRISHNAN, S. LANTERI, N. OLIVARES, AND R. PERRUSSEL 


we find that the worst values of r are along the imaginary axis. Finally, in Figures lg 
and lh, we see the effects of multiplying these real values of kh by 1 + i. The region 
of the complex plane where bad values of r are found changes significantly when kh is 
complex. 


Results on unisolvent stabilization 

We now turn to the question of how we can choose a value for the stabilization 
parameter r that will guarantee that the local matrices are not singular. The answer, 
given by a condition on t, surprisingly also guarantees that the global condensed HDG 
matrix is nonsingular. These results are based on a tenuous stability inherited from 
the fact nonzero polynomials are never waves, stated precisely in the ensuing lemma. 
Then we give the condition on r that guarantees unisolvency, and before concluding the 
section, present some caveats on relying solely on this tenuous stability. 

As is standard in all HDG methods, the unique solvability of the element problem 
allows the formulation of a condensed global problem that involves only the interface 
unknowns. We introduce the following notation to describe the condensed systems. 
First, for Maxwell’s equations, for any rj e Nh, let (E v , H r <) e Yj, x W denote the fields 
such that, for each K e 7*,, the pair (E ri \ K , H v \k) satisfies the local problem (11) with 
data J]\dK- That is, 

ik(E v , v)k - (Vx H' 1 , v)k + {tE v xh, fix ti)qk = (tt/ x h, fix n)sK, (15a) 

-(E n ,Vxw) K -ik(H 11 ,w) I< = -( rj,nxw ) dK , (15b) 

for all v e Y(K), w e Y(K ). If all the sources in (10) vanish, then the condensed global 
problem, for E e Nh takes the form 

a(E,rj)=0, Vr) € N h , (16) 

where 

o(A, v) = E (H A x n, rj)dK- 

KeT h 

By following a standard procedure [4] we can express a(-, ■) explicitly as follows: 

a( A, rj) = ^ ( H a x h, rj)g K + (( H A - H A ) x n, rj) dK 

KzT h 

= E ik(H A , H v ) k - (Vx H a , E v ) k + (rh x (ft x (A - E A )), 7j) dK 
KeT h 

- E ik(H A , H v )k ~ ik(E A , E V ) K + {rh x (A - E A ), n x E^) dK 
KeT h 

- {rh x (A - E A ), h x p) dK 

- E ik(H A , H v ) k - tk(E A , E v ) k - r(h x (A - E A ), n x (r/ - E n )) dK . 

KiT h 

Here we have used the complex conjugate of (15b) with w = H A , along with the definition 
of H A , and then used (15a). 
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Similarly, for the Helmholtz equation, let e 14 x H4 denote the fields such 

that, for all K e 7/,, the functions (u v \ki ^Ik) solve the element problem (4) for given 
data 0 = 77 . If the sources in (2) vanish, then the condensed global problem for 0 e M/, is 
written as 

6 ( 0 , 77 ) = 0 , V 77 e M h , (17) 

where the form is found, as before, by the standard procedure: 

KKv) = E {u A -n,v)dK 

KeT h 

= E ik(u A ,u ,1 ) K -ik(<j) A ,<j) v ) K -T(A-<j) A ,'ri-(j) ri ) dK . 

KeTh 

The sesquilinear forms a(-,-) and 6(-,-) are used in the main result, which gives suffi¬ 
cient conditions for the solvability of the local problems (11), (4) and the global prob¬ 
lems (16), (17). 

Before proceeding to the main result, we give a simple lemma, which roughly speaking, 
says that nontrivial harmonic waves are not polynomials. 

Lemma 1. Let p ^ 0 be an integer, 0 * k e C, and D an open set. Then, there is no 
nontrivial E e (V P (D)) 3 satisfying 

Vx(Vxi7) - k 2 E = 0 

and there is no nontrivial 0 € V P (D) satisfying 

A 0 + k' 2 (p = 0. 

Proof. We use a contradiction argument. If E f 0, then we may assume without loss of 
generality that at least one of the components of E is a polynomial of degree exactly p. 
But this contradicts k 2 E = Vx(Vx E) because all components of Vx(Vx E) are polyno¬ 
mials of degree at most p- 2. Hence E = 0. An analogous argument can be used for the 
Helmholtz case as well. □ 

Theorem 1. Suppose 

Re(r) + 0, whenever Im(7) = 0, and (18a) 

Im(/c) Re(r) ^ 0, whenever lm(k) + 0. (18b) 

Then, in the Maxwell case, the local element problem (11) and the condensed global 
problem (16) are both unisolvent. Under the same condition, in the Helmholtz case, the 
local element problem (4) and the condensed global problem (17) are also unisolvent. 

Proof. We first prove the theorem for the local problem for Maxwell’s equations. As¬ 
sume (18) holds and set E = 0 in the local problem (11). Unisolvency will follow by 
showing that E and H must equal 0. Choosing v = E, and w = H, then subtract¬ 
ing (lib) from (11a), we get 

(\\E\\k + \\H\\ 2 k) + 2ilm(E,WxH) K + t\\E x nf dK = 0, 
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whose real part is 

-ImO) (\\E\\ 2 k + \\H\\ 2 k ) + Re(r)||U x n\\ 2 dK = 0. 

Under condition (18b), we immediately have that the helds E and H are zero on K. 
Otherwise, (18a) implies E x n\g K = 0, and then (11) gives 

zkE - Vx H - 0, 

ikH + Vx E - 0, 

implying 

Vx(Vxh) = k 2 E. 

By Lemma 1 this equation has no nontrivial solutions in the space Y(K). Thus, the 
element problem for Maxwell’s equations is unisolvent. 

We prove that the global problem for Maxwell’s equations is unisolvent by showing 
that E - 0 is the unique solution of equation (16). This is done in a manner almost 
identical to what was done above for the local problem: First, set r) = E in equation (16) 
and take the real part to get 

E lm(*0 (II# lit + lldlt) - Re(r)||ft * (E - E)f m = 0. (19) 

!:--n 

This immediately shows that if condition (18b) holds, then the helds E and E[ are zero on 
c M 3 and the proof is finished. In the case of condition (18a), we have nx (E-E\qk) = 0 
for all K. Using equations (10), this yields 

[Vx(vxE)]\ K = k 2 E\ K , 

so Lemma 1 proves that the helds on element interiors are zero, which in turn implies 
E - 0 also. Thus, the theorem holds for the Maxwell case. 

The proof for the Helmholtz case is entirely analogous. □ 

Note that even with Dirichlet boundary conditions and real k, the theorem asserts the 
existence of a unique solution for the Helmholtz equation. However, the exact Helmholtz 
problem (1) is well-known to be not uniquely solvable when k is set to one of an infinite 
sequence of real resonance values. The fact that the discrete system is uniquely solvable 
even when the exact system is not, suggests the presence of artificial dissipation in HDG 
methods. We will investigate this issue more thoroughly in the next section. 

However, we do not advocate relying on this discrete unisolvency near a resonance 
where the original boundary value problem is not uniquely solvable. The discrete ma¬ 
trix, although invertible, can be poorly conditioned near these resonances. Consider, 
for example, the Helmholtz equation on the unit square with Dirichlet boundary condi¬ 
tions. The first resonance occurs at k = In Figure 2, we plot the condition number 

Q' ma .v/o'min of the condensed HDG matrix for a range of wavenumbers near the resonance 
k = vrV2, using a small fixed mesh of mesh size h = 1/4, and a value of r = 1 that sat¬ 
isfies (18). We observe that although the condition number remains finite, as predicted 
by the theorem, it peaks near the resonance for both the p - 0 and the p = 1 cases. We 
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(a) Degree p - 0 (b) Degree p = 1 

Figure 2. Conditioning of the HDG matrix for the Helmholtz equation 
near the first resonance k = i\\pl » 4.44. 

also observe that a parameter setting of r = -% that does not satisfy the conditions of 
the theorem produce much larger condition numbers, e.g., the condition numbers that 
are orders of magnitude greater than 10 10 (off axis limits of Figure 2b) for k near the 
resonance were obtained for p = 1 and r = To summarize the caveat, even though 
the condition number is always bounded for values of r that satisfy (18), it may still be 
practically infeasible to find the HDG solution near a resonance. 

Results of dispersion analysis for real wavenumbers 

When the wavenumber k is complex, we have seen that it is important to choose the 
stabilization parameter r such that (18b) holds. We have also seen that when k is real, 
the stability obtained by (18a) is so tenuous that it is of negligible practical value. For 
real wavenumbers, it is safer to rely on stability of the (un-discretized) boundary value 
problem, rather than the stability obtained by a choice of r. 

The focus of this section is on real k and the Helmholtz equation (1). In this case, 
having already separated the issue of stability from the choice of r, we are now free 
to optimize the choice of r for other goals. By means of a dispersion analysis, we now 
proceed to show that some values of r are better than others for minimizing discrepancies 
in wavespeed. Since dispersion analyses are limited to the study of propagation of 
plane waves (that solve the Helmholtz equation), we will not explicitly consider the 
Maxwell HDG system in this section. However, since we have written the Helmholtz 
and Maxwell system consistently with respect to the stabilization parameter (see the 
transition from (3) to (9) via (7)), we anticipate our results for the 2D Helmholtz case 
to be useful for the Maxwell case also. 

The dispersion relation in the one-dimensional case. Consider the HDG method (2) 
in the lowest order (p = 0) case in one dimension (ID) - after appropriately interpreting 
the boundary terms in (2). We follow the techniques of [1] for performing a dispersion 
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analysis. Using a basis on a segment of size h in this order u\ = 1, <f>\ - 1, 

1, 02 = 1) the HDG element matrix takes the form M = \ ] w here 


n = 


M n = 


ikh 

0 


0 

-ikh - 2 r 


Mi 2 - 


-1 +1 


T 


T 


M 2 1 - M\ 2 


Af 2 2 - 


-r 0 
0 -r 


The Schur complement for the two endpoint basis functions {0i,0 2 } is then 


S = 


1 

ikh 

1 

L ikh 




ikh + 2 r 




+ r 


1 

ikh 

1 

ikh 


t “ 


ikh + 2 r 

-2 

+ r 




ikh + 2 r ikh ikh + 2 r 

Applying this matrix on an infinite uniform grid (of elements of size /i), we obtain the 
stencil at an arbitrary point. If fij denotes the solution (trace) value at the 71 h point 
(j € Z), then the jth equation reads 


2% FIT 


T 


+ r ) + (0i-i + 0i+i)(' 


1 

ikh 


t “ 


= 0 . 


fikh ikh + 2r ) ' \ ikh ikh + 2r, 

In a dispersion analysis, we are interested in how this equation propagates plane waves 
on the infinite uniform grid. Hence, substituting tfj = exp (ik h jh), we get the following 
dispersion relation for the unknown discrete wavenumber k h : 


cos (k h h) ^ 




ikh ikh + 2 r 


)■( 


1 

ikh 




ikh + 2 r 


+ r 


Simplifying, 


k h h = cos 


-1 


1 - 


(kh.y 


( 20 ) 


2 + ikh(r + r _1 ) ) 

This is the dispersion relation for the HDG method in the lowest order case in one 
dimension. Even when r and k are real, the argument of the arccosine is not. Hence 

lm(k h ) + 0 , ( 21 ) 


in general, indicating the presence of artificial dissipation in HDG methods. Note how¬ 
ever that if r is purely imaginary and kh is sufficiently small, ( 20 ) implies that lm(k h ) = 0 . 

Let us now study the case of small kh (i.e., large number of elements per wavelength). 
As kh -* 0, using the approximation cos _1 (l - x 2 12) « x + x 3 /24 + ••• valid for small x, 
and simplifying ( 20 ), we obtain 

k h h -kh = - ^ ^ \ kh ) 2 + 0((kh) 3 ). ( 22 ) 

Comparing this with the discrete dispersion relation of the standard finite element 
method in one space dimension (see [ 1 ]), namely k h h - kh « 0((kh.) 3 ), we hnd that 
wavespeed discrepancies from the HDG method can be larger depending on the value of 
r. In particular, we conclude that if we choose t - ±i, then the error k h h - kh in both 
methods are of the same order 0((kh) 3 ). 
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Before concluding this discussion of the one-dimensional case, we note an alternate 
form of the dispersion relation suitable for comparison with later formulas. Using the 
half-angle formula, equation (20) can be rewritten as 

.M-l®!_I_ 

V 2 / \ikh(r 2 + 1) + 2r 

where c - cos(k h h/2). 



Lowest order two-dimensional case. In the 2D case, we use an infinite grid of square 
elements of side length h. The ffDG element matrix associated to the lowest order (p = 0) 
case of (2) is now larger, but the Schur complement obtained after condensing out all 
interior degrees of freedom is only 4x4 because there is one degree of freedom per 
edge. Note that horizontal and vertical edges represent two distinct types of degrees of 
freedom, as shown in Figures 5a and 5b. Hence there are two types of stencils. 

For conducting dispersion analysis with multiple stencils, we follow the approach in [7] 
(described more generally in the next subsection). Accordingly, let C\ and C2 denote 
the infinite set of stencil centers for the two types of stencils present in our case. Then, 
we get an infinite system of equations for the unknown solution (numerical trace) values 
-01^ and 0 2,p 2 at all p\ e C\ and p 2 e C 2 , respectively. We are interested in how this 
infinite system propagates plane wave solutions in every angle 9. Therefore, with the 
ansatz p^p :i = exp(i7y ■ pj) for constants a 3 (j = 1 and 2), where the discrete wave 
vector is given by 


K h = k h 


cos 9 
sin 9 


we proceed to find the relation between the discrete wavenumber k h and the exact 
wavenumber k. 

Substituting the ansatz into the infinite system of equations and simplifying, we obtain 
a 2 x 2 system F[ ] = 0 where 


F = 

and, for j = 1,2, 


2khr 2 ciC2 d\ (4r + ikh) + 2 khr 2 Ci 2 

d 2 ( 4 t + ikh) + 2khr 2 c 2 2 2 khT 2 c\C 2 


Cj = cos | -1 


,^-hkj^ , dj = 2 i(l - Cj) - rkh, ki=k h cos9, k, 2 =k h sm 6 . 
Hence the 2D dispersion relation relating k h to k in the HDG method is 


(24) 


det(F) = 0. 


(25) 


To formally compare this to the ID dispersion relation, consider these two sufficient 
conditions for det(F) = 0 to hold: 

2 (/c/?.) 2 r 2 c| + dj (2 rkh + i{kjh ) 2 ) = 0, 


for j = 1,2, 


(26) 
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where Aq = k cos 9 and k 2 = k sin#. (Indeed, multiplying (26)j by d 3 +\ (j = 1) or d 3 _\ 
(j = 2) and summing over j = 1,2, one obtains a multiple of det(F).) The equations 
in (26) can be simplified to 

2 i (kjh) 2 ( khr 

Cj 2 1 \(kjh) 2 + (kh) 2 t 2 -2ikhr 


3 = 1,2, 


(27) 


which are relations that have a form similar to the ID relation (23). Hence we use 
asymptotic expansions of arccosine for small kh , similar to the ones used in the ID case, 
to obtain an expansion for kb, for j = 1 , 2 . 

The final step in the calculation is the use of the simple identity 

/ \V2 


k h = (Aq fe ) 2 + (fc £) 2 


(28) 


Simplifying the above-mentioned expansions for each term on the right hand side above, 
we obtain 


k h h - kh - ‘V°s(49)+3 + m) + 3 

16 T 


(29) 


as kh -»• 0. Thus, we conclude that if we want dispersion errors to be 0((kh) 3 ), then we 
must choose 

r = ±-£\/cos(4d) + 3, (30) 


a prescription that is not very useful in practice because it depends on the propagation 
angle 6 . However, we can obtain a more practically useful condition by setting r to be 
the constant value that best approximates ±^iycos(4d) + 3 for all 0 ^ 6 ^ 7 t/2 , namely 



(31) 


These values of r asymptotically minimize errors in discrete wavenumber over all angles 
for the lowest order 2D HDG method. Note that for any purely imaginary r, (27) implies 
that kJf is real if kh is sufficiently small, so 


Im (k h ) = 0 , 


(32) 


thus eliminating artificial dissipation. 

We now report results of numerical computation of k h - k h (9 ) by directly applying 
a nonlinear solver to the 2D dispersion relation (25) (for a set of propagation angles 
9). The obtained values of the real part R ek h (9) are plotted in Figure 3a, for a few 
fixed values of r. The discrepancy between the exact and discrete curves quantifies the 
difference in the wave speeds for the computed and the exact wave. Next, analyzing 
the computed k h (9 ) for values of r on a uniform grid in the complex plane, we found 
that the values of r that minimize |kh - k h (9)h\ are purely imaginary. As shown in 
Figure 4, these r -values approach the asymptotic values determined analytically in (30). 
A second validation of our analysis is performed by considering the maximum error over 
all 9 for each value of r and then determining the practically optimal value of r. The 
results, given in Table 2, show that the optimal r values do approach the analytically 
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Exact wave speed 


Exact wave speed 

Discrete wave speed, x=1 


Discrete wave speed, x=1 

Discrete wave speed, x=-i 


Discrete wave speed, x=-i 

Discrete wave speed, t=0.1 +0.8i 


Discrete wave speed, x=0.1+0.9i 




Figure 3. Numerical wave speed R e(k h (6)) as a function of 9 for various 
choices of r. Here, k = 1 and h - 7 t/4 . 


kh 

Optimal r, 
Im(r) > 0 

Optimal r, 
Iui(t) < 0 

7 r/4 

o 

bo 

o 

> 

-0.931? 

7t/8 

0.837? 

-0.898? 

7t/16 

0.851? 

-0.882? 

vr/32 

0.859? 

-0.874? 

7r/64 

0.863? 

-0.871? 

vr/128 

0.865? 

-0.868? 

vr/256 

0 .866? 

-0.867? 


Table 2. Numerically found values of r that minimize \kh- k h \9)h\ for 
all 6 in the p - 0 case. 


determined value (see (31)) of ±?^ » ±0.866?. Further numerical results for the p - 0 
case are presented together with a higher order case in the next subsection. 

Higher order case. To go beyond the p - 0 case, we extend a technique of [7] (as 
in [9]). Using a higher order HDG stencil, we want to obtain an analogue of (25), which 
can be numerically solved for the discrete wavenumber k h = k h \9). The accompanying 
dispersive, dissipative, and total errors are defined respectively by 

e disp = max| R e(k h (9)) - k\, e dissip = max | Im(k h (0))|, 

, (33) 

^total = max | k h (9) - k\. 

Again, we consider an infinite lattice of h x h square elements with the ansatz that 
the HDG degrees of freedom interpolate a plane wave traveling in the 6 direction with 
wavenumber k h . The lowest order and next higher order HDG stencils are compared 
in Figure 5. Note that the figure only shows the interactions of the degrees of freedom 
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Angle 0 



0 7t/4 ti/2 

Angle 0 


(a) Im(T) > 0 


(b) Im(r) < 0 


FIGURE 4. The values of r that locally minimize |kh - k h h\ are purely 
imaginary. Here, (Im(r)) 2 is compared with asymptotic values (dashed 
lines). 




(d) 

FIGURE 5. Stencils corresponding to the shaded node types. (A)-(B): 
Two node types of the lowest order (p - 0) method; (C)-(F): Four node 
types of the first order (p = 1) method. 



corresponding to the 0 variable—the only degrees of freedom involved after elimination 
of the u and 0 degrees of freedom via static condensation. The lowest order method has 
two node types (shown in Figures 5a-5b), while the first order method has four node 
types (shown in Figures 5c-5f). For a method with S distinct node types, denote the 
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solution value at a node of the s th type, 1 ^ s ^ S, located at Ih e M 2 , by ifj s p With our 
ansatz that these solution values interpolate a plane wave, we have 


for some constants a s . 

Now, to develop notation to express each stencil’s equation, we fix a stencil within the 
lattice. Suppose that it corresponds to a node of the t th type, 1 ^ t ^ S, that is located 
at jh. For 1 ^ s ^ S', define J t)S = {h l 2 : a node of type s is located at (j + l)h} and, 
for l e J i)S , denote the stencil coefficient of the node at location (j + l)h by D ts j. The 
stencil coefficient is the linear combination of the condensed local matrix entries that 
would likewise appear in the global matrix of equation (17). Both it and the set J t)S are 
translation invariant, i.e., independent of j. Since plane waves are exact solutions to the 
Helmholtz equation with zero sources, the stencil’s equation is 

E E D t,s,l Vby+i = o. 

s=l 

Finally, we remove all dependence on j in this equation by dividing by e lk '^ h , so there 
are S equations in total, with the t th equation given by 

= » (34) 

Defining the S x S matrix F(k h ) by 

[F(k h )] Us = £ D t s l e ^ h [cosO, S inO]-lh^ 
lzJt,s 

we observe that non-trivial coefficients {a<j} exist if and only if k h is such that 

det(F(k h )) = 0. (35) 

This is the equation that we solve to determine k h for a given 6 for any order (cf. (25)). 

Results of the dispersion analysis are shown in Figures 3 and 6. These figures com¬ 
bine the results from previously discussed p - 0 case and the p - 1 cases to facilitate 
comparison. Here, we set k = 1 and h = 7 t/4 , i.e., 8 elements per wavelength. Figure 6 
shows the dispersive, dissipative, and total errors for various values of r e C. For both 
the lowest order and first order cases, although the dispersive error is minimized at a 
value of r having nonzero real part, the total error is minimized at a purely imaginary 
value of r. This is attributed to the small dissipative errors for such r. Specifically, the 
total error is minimized when r = 0.87? in the p = 1 case. This is close to the optimal 
value of r found (both analytically and numerically) for p = 0. This value of r reduces 
the total wavenumber error by 90% in the p - 1 case, relative to the total error when 
using t = 1. 



Im(i) lm(x) im(x) 
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Re(x) 


(a) Dispersive error, p = 0 



Re(x) 

(b) Dispersive error, p = 1 



Re(x) 

(e) Total error, p = 0 


Re(x) 

(f) Total error, p = 1 



Figure 6 . Normalized dispersive error edi sp /ed isp , dissipative error 
edissip/^dissip’ and total error e t.otai/e t ' 0 tai for various r e C. Here, k = 1, 
h = 7r/4, and , e dissi P an d e total denote the errors when r = 1, respec¬ 
tively. 
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Comparison with dispersion relation for the Hybrid Raviart-Thomas method. 

The HRT (Hybrid Raviart-Thomas) method is a classical mixed method [2, 3, 15] which 
has a similar stencil pattern, but uses different spaces. Namely, the HRT method for 
the Helmholtz equation is defined by exactly the same equations as (2) but with these 
choices of spaces on square elements: V(K) = Q p+ 1 jP ( K) x Q PiP+ 1 ( K ), W(K) = Qp(K), 
and M(F ) = V P (F). Here Qi, m (K ) denotes the space of polynomials which are of degree 
at most l in the first coordinate and of degree at most m in the second coordinate. The 
general method of dispersion analysis described in the previous subsection can be applied 
for the HRT method. We proceed to describe our new findings, which in the lowest order 
case includes an exact dispersion relation for the HRT method. 

In the p = 0 case, after statically condensing the element matrices and following the 
procedure leading to (25), we find that the discrete wavenumber k h for the HRT method 
satisfies the 2D dispersion relation 

( c i + c !) (2 (hk) 2 - 12) + c\c\ (4 (hk) 2 + 48) + ( hk) 2 - 24 = 0, (36) 

where Cj, as defined in (24), depends on k 1 -', which in turn depends on k h . Similar to the 
HDG case, we now observe that the two equations 


(2 (hkj) 2 + 12) c 2 + (hkj ) 2 - 12 - 0, j = 1,2, (37) 


are sufficient conditions for (36) to hold. Indeed, if lj is the left hand side above, then 
h(2c 2 + 1) + l 2 (2c\ + 1) equals the left hand side of (36). The equations of (37) can 
immediately be solved: 


hkj = 2 cos 1 


12 -{hkjf 


1/2 


2 (hkj) 2 + 12) 

Hence, using (28) and simplifying using the same type of asymptotic expansions as the 
ones we previously used, we obtain 


k h h -kh = - COS ^ + 3 (kh) 3 + 0((khf) (38) 

as kh -> 0. Comparing with (29), we find that in the lowest order case, the HRT method 
has an error in wavenumber that is asymptotically one order smaller than the HDG 
method for any propagation angle, irrespective of the value of r. 

To conclude this discussion, we report the results from numerically solving the nonlin¬ 
ear solution (36) for k h (9) for an equidistributed set of propagation angles 9. We have 
also calculated the analogue of (36) for the p = 1 case (following the procedure described 
in the previous subsection). Recall the dispersive, dissipative, and total errors in the 
wavenumbers, as defined in (33). After scaling by the mesh size h, these errors for both 
the HDG and the HRT methods are graphed in Figure 7 for p - 0 and p = 1. We find 
that the dispersive errors decrease at the same order for the HRT method and the HDG 
method with r = 1. While (38) suggests that the dissipative errors for the HRT method 
should be of higher order, our numerical results found them to be zero (up to machine 
accuracy). The dissipative errors also quickly fell to machine zero for the HDG method 
with the previously discussed “best” value of r = b/3/2, as seen from Figure 7. 



|kh—k h| |lm(k h)| |kh-Re(k h h)| 


22 


J. GOPALAKRISHNAN, S. LANTERI, N. OLIVARES, AND R. PERRUSSEL 



(c) Dissipative error, p = 0 (d) Dissipative error, p = 1 




(e) Total error, p = 0 (f) Total error, p = 1 


-®- HDG with t = 1 

-Order (kh) 2 reference 

-■- HDG with x = 0.866 i 

Order (kh) 3 reference 

-*- HRT 

.Order (kh) 4 reference 


Order (kh) 5 reference 


Figure 7. Convergence rates as kh -*■ 0 
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Conclusions 

These are the findings in this paper: 

(1) There are values of stabilization parameters r that will cause the HDG method 
to fail in time-harmonic electromagnetic and acoustic simulations using complex 
wavenumbers. (See equation (5) et seq.) 

(2) If the wavenumber k is complex, then choosing r so that Re(r)lm(A;) ^ 0 guar¬ 
antees that the HDG method is uniquely solvable. (See Theorem 1.) 

(3) If the wavenumber k is real, then even when the exact wave problem is not well- 
posed (such as at a resonance), the HDG method remains uniquely solvable when 
Re(r) + 0. However, in such cases, we found the discrete stability to be tenuous. 
(See Figure 2 and accompanying discussion.) 

(4) For real wavenumbers k , we found that the HDG method introduces small 
amounts of artificial dissipation (see equation (21)) in general. However, when r 
is purely imaginary and kh is sufficiently small, artificial dissipation is eliminated 
(see equation (32)). In ID, the optimal values of r that asymptotically minimize 
the total error in the wavenumber (that quantifies dissipative and dispersive er¬ 
rors together) are r = ±i (see equation (22)). 

(5) In 2D, for real wavenumbers k, the best values of r are dependent on the propa¬ 
gation angle. Overall, values of r that asymptotically minimize the error in the 
discrete wavenumber (considering all angles) is r = ±b/3/2 (per equation (31)). 
While dispersive errors dominate the total error for r = b/3/2, dissipative errors 
dominate when r = 1 (see Figure 7). 

(6) The HRT method, in both the numerical results and the theoretical asymptotic 
expansions, gave a total error in the discrete wavenumber that is asymptotically 
one order smaller than the HDG method. (See (38) and Figure 7.) 
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